##
##The following code replicates Figure A.3 from the supplemental appendix to our main paper
##

#clear memory  
rm(list=ls())

#load package  
library(TeachingDemos)
library(foreign)
library (gplots)

#read.ziop data
data.ziop<-read.dta("Global.RMSEs.ZIOPDGP.dta")

#set up graphing parameters
par(cex.lab=1.1)
par(cex.axis=1.1)
par(cex.main=1.1)
par(mfrow=c(1,2))

########Beta 1 Graph##########

par(mar=c(5,4,2,2)+0.1)

ziop.dgp<-cbind(data.ziop$RMSE_mean_OP,data.ziop$RMSE_mean_ZIOP,data.ziop$RMSE_mean_ZIOPC)
colnames(ziop.dgp) <- c("OP", "ZIOP","ZIOPC")
ziop.dgp<-as.data.frame(ziop.dgp)

#now graph:
plot_colors <- c("red","green","blue")

#plot OP
plot(ziop.dgp$OP, type="o", pch=21, col=plot_colors[1], ylim=c(0,.34), axes=FALSE, ann=FALSE,cex = 1.5)

#create axis
axis(1, at=1:7, lab=c("0%","10%","30%","50%","70%","90%","100%"))
axis(2, at=c(0.0,0.1,0.2,0.3,0.4,0.5),lab=c("0.0","0.1","0.2","0.3","0.4","0.5"))

# Create box around plot
box()

# Graph MiOP
lines(ziop.dgp$ZIOP, type="o", pch=22, lty=2, col=plot_colors[2],cex = 1.5)

# Graph MiOPC
lines(ziop.dgp$ZIOPC, type="o", pch=24, lty=3, col=plot_colors[3],cex = 1.5)

#label axis
title(xlab= "Proportion of Non-Inflated Observations")
title(ylab= "Average RMSE")
title(main= "Marginal Effect Estimates (ZiOP DGP)")

legend(1.2, .33, legend=c("OP","ZIOP","ZIOPC"), cex=1, col=plot_colors, 
   pch=c(21,22,24), lty=1:3)


###############################################################################################

#read.ziopc data
data.ziopc<-read.dta("Global.RMSEs.ZIOPCDGP.dta")

#set up graphing parameters
par(cex.lab=1.1)
par(cex.axis=1.1)
par(cex.main=1.1)

########Beta 1 Graph##########

par(mar=c(5,4,2,2)+0.1)


ziopc.dgp<-cbind(data.ziopc$RMSE_mean_OP,data.ziopc$RMSE_mean_ZIOP,data.ziopc$RMSE_mean_ZIOPC)
colnames(ziopc.dgp) <- c("OP", "ZIOP","ZIOPC")
ziopc.dgp<-as.data.frame(ziopc.dgp)

#now graph:
plot_colors <- c("red","green","blue")

#plot OP
plot(ziopc.dgp$OP, type="o", pch=21, col=plot_colors[1], ylim=c(0,.34), axes=FALSE, ann=FALSE,cex = 1.5)

#create axis
axis(1, at=1:7, lab=c("0%","10%","30%","50%","70%","90%","100%"))
axis(2, at=c(0.0,0.1,0.2,0.3,0.4,0.5),lab=c("0.0","0.1","0.2","0.3","0.4","0.5"))

# Create box around plot
box()

# Graph MiOP
lines(ziopc.dgp$ZIOP, type="o", pch=22, lty=2, col=plot_colors[2],cex = 1.5)

# Graph MiOPC
lines(ziopc.dgp$ZIOPC, type="o", pch=24, lty=3, col=plot_colors[3],cex = 1.5)

#label axis
title(xlab= "Proportion of Non-Inflated Observations")
title(ylab= "Average RMSE")
title(main= "Marginal Effect Estimates (ZiOPC DGP)")

legend(1.2, .33, legend=c("OP","ZIOP","ZIOPC"), cex=1, col=plot_colors, 
   pch=c(21,22,24), lty=1:3)

